Ferromagnetic phase transition in a Heisenberg fluid: Monte 
Carlo simulations and Fisher corrections to scaling 



I. M. Mryglod, 1,2 , I. P. Omelyan, 1 and R. Folk 2 
1 Institute for Condensed Matter Physics, 1 Svientsitskii Str., UA-79011 Lviv, Ukraine 
2 Institute for Theoretical Physics, University of Linz, A-4040 Linz, Austria 

(February 1, 2008) 

Abstract 



The magnetic phase transition in a Heisenberg fluid is studied by means of 
the finite size scaling (FSS) technique. We find that even for larger systems, 
considered in an ensemble with fixed density, the critical exponents show de- 
viations from the expected lattice values similar to those obtained previously. 
This puzzle is clarified by proving the importance of the leading correction 
to the scaling that appears due to Fisher renormalization with the critical 
exponent equal to the absolute value of the specific heat exponent a. The 
appearance of such new corretions to scaling is a general feature of systems 
with constraints. 
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Monte Carlo (MC) simulations of finite systems near phase transitions have received 
considerable attention in recent years Of notable current interest are continuum spin 
fluid models [Q] which are considered as a first step towards the modelling of ferrofluids |3| 
and adsorption surface phenomena [|J . Several important results, which have both theoret- 
ical and experimental interest, were obtained for such models. For example, it was found 
that the phase diagrams for spin fluids due to the interplay between spin and translational 
degrees of freedom are much more complicated || than in non-magnetic liquids. Magnetic 
ordered phases can exist both in gas and liquid states. By applying an external magnetic 
field, one can shift significantly the locus of the gas-liquid transition || and change the dy- 
namic properties ||; both the static and dynamic properties in the models discussed show 
differences from the non magnetic fluid and the magnetic lattice model. 

One important question in this context is whether the magnetic transition in a Heisen- 
berg fluid belongs to the same universality class as the corresponding transition in the lattice 
model. On general grounds (annealed systems) J/J, the lattice universality class is expected. 
In Ref. Q using MC method, a novel set of critical exponents was found that were in dis- 
agreement with the expected results. Similar disagreements were later obtained for two- and 
three-dimensional (2d and 3d) Ising fluids , where Fisher renormalized exponents were 
expected 0. In all the cases mentioned, a weak dependence of universal quantities on the 
density of particles n = N/V and systematic deviation from the predicted critical exponents 
were observed in the MC simulations. The general conclusion was that computer simulations 
were strongly affected by nonlinear crossover effects, which hide the true asymptotic critical 
behavior, giving only effective critical exponents. 

The goal of the present study is to resolve this puzzle by performing a new series of MC 
simulations for a Heisenberg fluid, considering larger finite size systems, and to compare the 
results obtained with the previous data || . Contrary to the belief that the true asymptotic 
critical behavior would be observed, deviations from the expected exponents remained. It 
is the aim of this Letter to show that a new correction term ]7[ (Fisher correction term) has 
to be taken into account in the FSS even for very large systems. 
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Let us consider a classical system, composed of N magnetic particles of mass m and 
described by the Hamiltonian P,R| 



N IN 

mv. 



# = E ^ + E [ $ (^) - Jfa) Bi'Bj] , (1) 



i=l i<j 

where rj and Vj denote the position and velocity of particle % carrying spin Sj. In our 
MC study, the liquid subsystem potential $(ry) was chosen to be of soft-core-like form, 
$(r) = Ae[(a/r) 12 — (cr/r) 6 ] + e at r < 2 1 / 6 cr and $(r) = otherwise, and the exchange 
integral J{rij) > 0, describing spin interactions, was modelled by the Yukawa function, 
J(r) = e(cr/r) exp[(a — r)/<r]. The function J(r) was truncated at R = 2.5a and shifted 
to zero at the truncation point (this avoids force singularities during MD calculations [|TT|j). 
Staying within the classical approach we consider Sj as a three-component continuous vector 
with a fixed length |sj| = 1. 

The simulations were carried out in the basic cubic box V = L 3 (employing periodic 
boundary conditions) at the reduced density of n* = Na 3 /V = 0.6 for a reduced core in- 
tensity of e/e — 1. The number of particles N were taken as iV =125, 256, 512, 1000, 2048, 
4000, 8000, and 16384. The simulations have been performed for five values of reduced tem- 
perature T*, T* = k B T/e =2.000, 2.025, 2.050, 2.075, and 2.100. The system was allowed to 
achieve equilibrium for 100 000 iV attempted moves. The total number of trial moves per 
particle (cycles) performed in the equilibrium state was 1000 000. The canonical averaging 
over the system was carried out using a biasing scheme [12|,[13| for sampling orientational de- 



grees of freedom. To minimize computational costs, the cell list technique |13j was employed 
in handling the interparticle interactions. 

The critical properties of a system in the thermodynamic limit may be extracted from 
the behavior of finite size systems by examining the size dependence of thermodynamic 
quantities p|,[14||. According to the FSS theory, various thermodynamic quantities can be 



written in a scaling form q(L,T) = L Xq l v Q(z), where L is a linear length of system 



is a critical exponent of the quantity q, and z = tL 1 ^ is the temperature scaling variable 
with t = (T — T c )/T c (T c and v are the bulk critical temperature and critical exponent of 



correlation length £). Because we are interested in zero-field properties, only the scaling 
variable z appears in the scaling function Q. 

There are several methods to determine the critical temperature T c . One of them is the 
Binder crossing technique |I| formulated for the fourth-order cumulant U4 = 1 — M4/(3M|), 
where Mi = (m l ), m = |m|, and m = iV _1 J2iLi s i- This method does not need any 
assumptions about critical exponents, but for small systems the position of intersection 
points between any two curves U4, related to systems with lengths L and L', depends 
usually on L and L', because of corrections to FSS. We have estimated the value of T c as 
the average over the cross-section temperatures T cross (L, L'), found for systems with L = Li 
and V = L i+1 (i > 3), giving T c = 2.055 ± 0.001. Here and below the increasing subscript 
% denotes increasing numbers of particles iVj from the set considered. In a similar 

way, we found Ul = 0.618 ± 0.003. The same estimate for T c has been found by using the 
crossing technique for the function £(Li,T)/Li (see, e.g., ||15|| ) within the phenomenological 
renormalization group scheme. Using a more precise method for extracting T c [jXJ] from the 
values T cross (L, V), obtained for the Binder parameter U4 with the fixed ratio L4/L1 = 
L5/L2 = L 7 /L 4 = L 8 /L 5 = 2, we have obtained T c = 2.054 ± 0.001 and Ul = 0.619 ± 0.002. 

An alternative method, proposed in [|16||, allows one to estimate simultaneously both the 



critical temperature T c and the critical exponent v within the same series of calculations. 
The main idea of this approach (the scanning technique) is to look for a quantity-independent 
slope of the set of functions VJ with I = 1, 2, . . . , 6, all of which have similar scaling behavior. 
These functions are defined via the derivatives Ki = d Mi/ 0(3 (see, for details, Ref. [16|]). The 
results are shown in Fig. 1, so that we have got T c = 2.057 ±0.001 and 1/z/ = 1.396 ±0.006. 
Note that within the scanning technique, no corrections to scaling have been taken into 
account. 

We have also used other known methods [|1] to estimate T c . One of these (the shifting 
technique) is based on the analysis of the size- dependent shift of a peak T pea k(L), observed 
in some thermodynamic quantities (e.g., specific heat CV, susceptibility x, derivatives K L ). 
If corrections to scaling are neglected, the location of the peak T pea k(L) has the general 



form T pca k(L) = T c + AL^ 1 ^, where A is a quantity- dependent constant. Note that in 
order to determine T c one has to estimate accurately the exponent v as well as the values 
7peak(£)- In the temperature range considered, well-defined peaks for all the sizes L, have 
been observed for the functions dM/d(3, U 3 = (M 3 - 3M 2 Mi + Mf)/[M 1 (M 2 - M 2 )], and 
Xs = N 2 {M 3 -3M 2 M 1 + Mf) jrSyni- For all these cases, our estimate of T c ~ 2.058 ±0.002, 
found with 1/V = 1.396 for the five largest sizes L, is in agreement with the scanning 
technique but not with the crossing technique. Moreover, the dependence T peak (L) versus 
L v l v showed a pronounced curvature for smaller system sizes L. Hence, the first puzzle 
uncovered in our study is the disparity in the estimates for T c found using two types of 
standard FSS techniques, namely, (i) the crossing technique for Binder parameter as well 
as for the correlation length, and (ii) the scanning and shifting techniques. This disparity 
could not be entirely explained by the error bars and indicates strong crossover effects. Note 
that corrections to scaling were completely neglected in the methods of type (ii). In order to 
clarify the cause of the difference found for T c , the scanning technique was used for the four 
largest systems only, giving T c = 2.054 ± 0.001 in agreement with the result of the crossing 
method. However, despite the improvement in T c , the value 1/V = 1.312 ± 0.007 was a 
deterioration. This is another puzzle which needs explanation. Taking everything together 
these results support the presence of strong crossover effects in the system considered. 

Knowing T c , one can then estimate the critical exponents again using the FSS theory 
KJ. We have calculated the exponent ratios (3/v and 7/1/ the FSS behavior of M(L, T c ) and 
the magnetic susceptibility x(L,T c ). For T c = 2.057, we found (3/v — 0.544 ± 0.015 and 
j/iy — 1.90 ± 0.03, respectively. Nearly the same estimates were obtained using other FSS 
methods. Taking T c = 2.054 for the four largest sizes L i; we obtained (3/v = 0.520±0.008 and 
7/1/ = 1.87±0.03. The results for both choices T c = 2.057 and T c = 2.054 are summarized in 
Table 1 in the first and second lines, respectively, beneath the Fluid 2 heading. The estimates 
for 1/V were found by the scanning technique. 

Comparing our results with the previous ones ||, we conclude that: (i) the ratios of 
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critical exponents (3/u and 7/V found in our study are closer to the values known for the 
lattice model Jl6yi8|1 ; (ii) the critical exponent v is extremely sensitive to the estimate of 
critical temperature T c used; (iii) even for larger systems, which this study considers, a 
systematic deviation from the lattice exponents is seen that cannot entirely be justified by 
the error bars; and (iv) the disparity in estimates found for the critical temperature T c , using 
the crossing and shifting techniques, has no explanation within the standard FSS approach. 
Hence, our data have to be considered as results for effective exponents and one can expect 
that the true asymptotic behavior would be visible only for much larger systems. If non- 
asymptotic crossover effects are considered, one may think of the presence of the Wegner 
correction term. However, this is expected to be negligible for our largest system sizes. 
One has to ask, therefore, what the reason is for such a strong crossover in the system we 
considered, compared to the lattice model. 

In order to investigate this problem in more detail let us recall an idea encountered in 
the Fisher renormalization J7[ for a system under thermodynamic constraint. According 
to this idea, the critical singularities in the grand canonical ensemble, with fixed chemical 
potential /i, may be different from those describing the system in the canonical ensemble 
with fixed density n. One has to performed the corresponding Legendre transformation 
carefully, taking into account the properties of singular functions in the grand canonical 
ensemble. In particular, this gives the well-known relation (see, e.g., Eq. (2.38a) in 0) 

T = a t x "{l + a 1 t A "), (2) 

which connects the reduced temperature scales t and r in the two different ensembles with 
fixed n and /i, respectively. The values of x a and A a in (0) depend on the sign of the 
specific heat critical exponent a, and are equal to (1 — a) -1 or 1, and a(l — a)~ l or —a for 
a positive or negative, respectively. It is seen already from (0) that independent of whether 
Fisher renormalization changes the critical exponents in the ensemble with fixed n, a new 
type of corrections to scaling appears in the canonical ensemble. These corrections, being 
proportional to a, must not be confused with Wegner corrections to scaling; and because of 
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the smallness of a in the Heisenberg universality class, they have to be taken into account 



within the FSS analysis. We note also that Eq. (0) is not the only source |H| for the 
appearance of new corrections (as was assumed, e.g., in PUfl). There is another reason, 
which also follows from thermodynamics. For example, using hyperscaling relations for the 



critical exponents, it can be easily proved [HJ that the second term in the known expression 
XTn = (dM/dh) Tn = {dM/dh)^-{dM/dn) T (dN/d/j,)' 1 (dN/dh)^ produces an additional 
correction to the magnetic susceptibility xr,n with an exponent proportional to a. Hence, 
these new corrections to scaling cannot be included in the standard FSS by means of simple 
rescaling of t, which follows from (H) and as was proposed in [2(| . 



In order to prove our predictions and to estimate the range of asymptotic behavior in 
which the new correction can be neglected, we have performed additional calculations. In 
Fig. 2 the results, obtained for the temperatures T pea ^(L), where the maximums of the 
functions 8M/dj3 and x(T, L) are located, are shown for different sizes of L. These results 
have been fitted (dashed lines) for the four largest system sizes by using the expression 

T peak (i:) = T C + AN- 1 / 3 " (l + SiV-H/ 3 ") , (3) 



with the values of 1/u and a known for the lattice model |T^,[TB|. The quantity-dependent 



constants A and B were then estimated. It is seen in Fig. 2 that: (i) the fitting curves 
are in rather good agreement with the MC data obtained even for smaller values of L; 
(ii) such a simple procedure allows one to understand the strong deviation from the linear 
dependence T pea k(L) = T c + AN~ 1 / 3u that follows from (|3|) when the correction to scaling is 
neglected; and (iii) the disparity in estimates found for T c within the crossing and shifting 
techniques can be explained. Using the fitting procedure, described above, we have found 
that the estimate T c = 2.055±0.001 gives a rather good fit for all the data T peak (L), obtained 
from the maximum positions of dM/d/3, x{T,L), tJ^(T,L), and x 3 (T, L). Another finding 
was that in contrast to the strong quantity-dependence of A, the parameter B in ([|) is 
almost independent of the quantity considered |^TJ. Note that if the rescaling relation (^) is 
considered as the unique reason for the appearance of the new correction, then the parameter 



B is quantity-independent. From the fitting procedure it has been found that B ~ 1.3 ± 0.2 
for all the five sets of T pea ^(L) studied. Having the value of B, we can then estimate the 
minimal number of particles iV m i n such that the correction term in (|3|) can be neglected if 
iV > iV m i n . This gives iV m i n ~ 10 8 (then the relative contribution of the second term in the 
bracket of (|3|) is less than 0.5), and, therefore, it is clear why the true asymptotic behavior 
could not be observed earlier ||, or in our MC study. Finite systems with iV > A^ min have 
not been considered so far in MC simulations. Hence, only the effective exponents could be 
studied for smaller size systems. 

In conclusion, we note that if the absolute specific heat exponent a is small enough, 
Fisher corrections to scaling discussed are very important in models with constraints. This 
holds at d = 3 for the Ising, the XY and the Heisenberg classes of universality, as well as 
for other systems. In particular, we are convinced that the problems found in the Ising fluid 



10] have the same origin. In this respect it is also worth referring the reader to Refs. [22 



where, within the e-expansion scheme, it was proven analytically that the leading correction 
to scaling in a compressible Heisenberg magnet as well as in a randomly diluted, weakly 
inhomogeneous Heisenberg model is equal to —a, and this supports our conclusions. More 
detailed results including the determination of the values of the asymptotic exponents will 
be given elsewhere. 
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work was supported by the Fonds zur Forderung der wissenschaftlichen Forschung under 
Project No. P12422-TPH. 
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TABLE 1. Summary of results with n being the reduced density. In the first three rows 
(denoted as Fluid 1 ) the results from || are gives. In the last row the universal quantities 



known for the lattice Heisenberg model [|l6| , |T8f are presented. Our results are shown in the 
rows denoted as Fluid 2 . 





u 4 


1/u 


P/v 


l/v 


Fluid 1 










n=0.4 


0.613 


1.35(5) 


0.55(2) 


1.86(3) 


n=0.6 


0.608 


1.41(3) 


0.56(2) 


1.85(1) 


n=0.7 


0.605 


1.42(3) 


0.55(2) 


1.84(3) 


Fluid 2 










n=0.6 


0.619 


1.40(1) 


0.54(2) 


1.90(3) 






1.31(1) 


0.52(1) 


1.87(3) 


Lattice 


0.622 


1.421(5) 


0.514(1) 


1.973(2) 
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FIGURE CAPTIONS 



FIG. 1. Quantity dependence of scanning results for the functions VJ possessing the same 
scaling properties |1£]]. The horizontal line for T c = 2.057 is drawn at \jv — 1.396. 



FIG. 2. Size dependence of the maximum locations in the derivative DM/ 8(3 (triangles) 
and susceptibility x(T,L) (diamonds). The results of fitting to MC data (dashed curves) 
are found using (§) with \ jv and a known for the lattice model. 
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